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ABSTRACT 

This research focuses on the computational dy- 
namics of flexible constrained multibody systems. At 
first a recursive mapping formulation of the kinemat- 
ical expressions in a minimum dimension as well as 
the matrix representation of the equations of mo- 
tion are presented. The method employs Kane’s 
equation, FEM and concepts of continuum mechan- 
ics. The generalized active forces are extended to 
include the effects of high temperature conditions, 
such as creep, thermal stress and elastic-plastic de- 
formation. The time variant constraint relations for 
rolling/contact conditions between two flexible bod- 
ies are also studied. The constraints for validation 
of MBS simulation of gear meshing contact using 
a modified Timoshenko beam theory are also pre- 
sented. The last part deals with minimization of 
vibration/deformation of the elastic .beam in multi- 
body systems making use of time variant boundary 
conditions. The above methodologies and computa- 
tional procedures developed are being implemented 
in a program called D YAM US. 


reference frame of body j to element i, respectively. 
N is the shape function matrix, p denote the nodal 
coordinates, and n a set of unit vector fixed in R ( 
see reference [l]-[2] for more detail). 

The velocity of element i of body j found by differ- 
entiation of the above equation can be expressed as 


v Jt = {i} T {V J ’}+{j} T {V/}+{C} T {V/}+{p} T {V/'} 

( 2 ) 

Four arrays are identified in the velocity expres- 
sion and found to take a special form. Note that 
x represent the rigid body rotation between adjacent 
bodies. The partial derivative of the element velocity 
yield the following 


[ ([^] + [5 <3 ])[5 10 ] 

([^] + [5 <a ])[5 20 ] 


\y jt ] = m 


([Sr, j + mp*] 


( 3 ) 


KINEMATICS OF FLEXIBLE TREE-LIKE 
SYSTEMS 


0 


An explicit matrix representation of the partial 
velocities and partial angular velocities for tree-like 
structures is given below. Consider a flexible body in 
a MBS discretized into P elements. Let the position 
vector to an arbitrary element i of body j w.r.t. a 
fixed reference frame R be given by 

H{ j) H(j) 

h - 0 h-0 

+ ({r,J T + {*dW)P i# ]H*} (1) 

where 5 denotes the shift matrix, g, ( and r represent 
the body vector, the translation vector between ad- 
jacent bodies, and the position vector from the local 


where W is a transformation matrix used to isolate 
the generalized coordinate derivatives from the gen- 
eralized speeds. S qx 5^, S r and S JX are skew matri- 
ces corresponding to g, C, r and P> respectively. The 
partial velocity array associated with element defor- 
mation is given by 


[vn = 


IN) T [S>°] 


0 


( 4 ) 
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{*»} = E E / 

3 * 


The partial velocity array associated with q is ex- 
pressed as 

m 

[5‘°] 


[V?] = 


[S 


i -i*°i 


( 5 ) 


L o J 

The bodies of the above block matrices could be 
achieved through a budgeting procedure where a 
master block is first developed then the rest of the 
arrays are formed through a partition and mapping 
technique see Table 1. 


+ EE/ (s) 

It is important to note at this stage how the kine- 
matical expression form the bulk of all computations. 
In the above equations m Jt denotes the mass of ele- 
ment i in body j, I }t the tensor dyadic, f }t the force 
vector array acting on element i of body j, M )X the 
corresponding moment array and b Jt the surface trac- 
tion contribution vector. 


EQUATIONS OF MOTION 

The governing equations of motion for flexible 
multibody systems can be expressed as 


[M}{y) + [<?]{»} + [K}{y] = {F} (6) 

where [M] denotes the generalized mass matrix com- 
posed of 9 submatrices of the form 

( 7 ) 

The generalized mass is symmetric and the compo- 
nents of M XJ come directly from the kinematic bank 
of partial velocities and angular velocities of ele- 
ments. The other mass components have similar ex- 
pressions. Similarly, we can write the dynamic damp- 
ing matrix, generalized stiffness matrix and force vec- 
tor in a partition form with its components expressed 
as 


CONSIDERATION OF HIGH TEMPERA- 
TURE, CREEP AND ELASTIC-PLASTIC 
DEFORMATIONS 

The modeling of time-dependent forces resulting 
from deformable bodies when subjected high tem- 
perature conditions can be of interest in many en- 
gineering applications, which include creep, thermal 
stress, thermal shock, etc.. Many researchers studied 
material nonlinearities, in which some problems are 
solved, other still remain to be issues of concern. 

The effects of temperature, creep and thermal 
stress and thermal shock can be included in the third 
term of the generalized force (see reference [2]) 


{*} = E E/ + [«*]{ 


3 i *3' 


EE/ w 


+ {fr} (10) 


> ■ •> 


where the last part {Ft} brings in the contribution 
from the effects of temperatures and material nonlin- 
earities 


= E E/ k 


]([/,. M 


<*V} = EE/ + 


3 * 


+ [n' i ][J,.][u>'] T )}<*® (8) 


and 


+*W)] + *M T [-DKT'}}<it> (11) 
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The nonlinearities including geometric nonlinearity 
and material nonlinearity can be considered in the 
stiffness matrix. So does the elastic- pi as tic deforma- 
tion. The material property matrix is given by 

[D]=[D e ] + [D p ] 

where [jD«] denotes the elastic part. The second part 
[Dp\ is the contribution from the plastic deformation 


While the time variant contact conditions can be 
considered as a set of constraints which can be holo- 
nomic or nonholonomic. Some constraint equations 
which do not contain prescribed motion terms can be 
factorized to minimize the dimension of the equations 
of the system. For the case of two flexible bodies with 
one rolling without slipping on the other as shown in 
Figure 1, we can write in R the following position 

vector^ 


[D P ] = [DJ- 




P„, = P»-1 +«n + r " 


(17) 


( 12 ) 


where 


TIME VARIANT BOUNDARY 

CONSTRAINT CONDITIONS 

For the time variant boundary conditions, finite 
difference method can be used to account for the rate 
of change of mode shape. Consider the modal trans- 
formation 

M = [*}{p] ( 13 ) 


<l n — T n -l ,c -f - 1 ,c — ( r nc + ttnc) (1®) 

Differentiation of equation (17) yields the constraint 
equations at the velocity level 

[J]{v] = M ( 19 ) 

where 


Differentiation of the above equation yields 


{n} = [*M + (*Kp} ( 14 ) 


{y} = [ * T C T p t ] T (20) 


When substituting the nodal displacement with the 
nodal coordinates and taking into consideration the 
effects of [4], then at t = U the new terms coming 
from the previous and newly computed mode shapes 
at t = t] are seen in [C] and [K] as f2),l3] 

[Cnz] = [<? 33 ] + 0^) 

and 

[Kns] = [Kn]+^[M»]([*i] 

- 2[$,_i] + [*i_a]) + £j[C»]([*i] - [*<-»]) ( 16 ) 

The method developed above has a wide range of 
applications for which one can easily see and analyze 
its dynamics. 


{ 9 } = [5 n_1 ■"][O n_1,n ]({r n } + [Nn]{pn*}) (21) 

[J] is a Jacobi matrix and a function of generalized 
coordinates and velocities. 

In the dynamics of MBS for the case when one 
flexible body is rolling on another, equations (19) and 
(6) extended with XJ T are solved together. The time 
history of the system allows us to systematically up- 
date the contact position and the reevaluation of the 
Jacobi matrix J. 

DYNAMICS OF GEAR MESHING TEETH 

For validation of the results obtained by multibody 
dynamics code which utilize FEM, a modified Timo- 
shenko beam theory is presented to analyze the dy- 
namics of gear meshing teeth in rotorcraft systems. 
The acting position, direction and magnitude of the 
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external forces are assumed to time variant. The 
meshing tooth is Considered as a cantilever beam, as 
shown in Figure 2, where the inertia force due to 
the large rotation of the tooth base, as well as the 
external equivalent axial force and moment are all 
included in the equation of motion.^ 




d 2 w 1 r , ,d 2 w , . 


+ U^W [eA{x) ^ + fix ' t) + P{x ’ t) iS ] 

+ e A(*)^+/(*,t) + J > (M)0 = 0 (22) 


For the assumed model, the boundary conditions are 
given by: 

At the fixed end z = 0, 


V>(0,t) = w(0,t) = 0 > (23) 


At the free end x — /, 

V{l,t)=kA{l)G{^^--i,(l,t)} = 0 (24) 


M(l,t) = [EI(l)^f^+m(l,t)] = 0 (25) 

A solution to the above proposed model will result 
in prediction of contact forces or dynamic loading on 
gear teeth. 

MINIMIZATION OF VIBRATION IN 
ELASTIC BEAMS 

The minimization of vibration (deformation) of 
flexible bodies in mechanical systems is a major con- 
cern in dynamics and control. What follows are pro- 
cedures used to minimize vibration in elastic beams. 
The elastic beam is modeled in two ways: one has a 
movable support not to exceed the lower tip, whereas 
the other treats the body as a hollow beam with a 
moving mass. 

Equation of motion for the model used to minimize 
vibration of the flexible beam, as shown in Figure 3, 
is given by' 2 ) 


Ar* v M2 


+ gcosO -f Cq -f xO) = 0 (26) 


Laplace transform gives 


4 

Y= Y^c.e riX + ^{L(x,s)+m[sMx) + /,(*)]} 

(27) 

The functional used to minimize vibration of the 
beam is 



F(z, zo , t)dt 


Euler-Lagrange equation 


(28) 


BF 

dxo 


d BF 
dt Bx o 


) = 0 


(29) 


is used to solve for the problem at hand. 

The solution for optimum positioning conditions 
is time variant and yields minimum deflection at the 
proposed location of the beam. 

* 
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